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Abstract 

Covariant vectors, Lyapunov vectors, or Oseledets vectors are increasingly being 
used for a variety of model analyses in areas such as partial differential equations, 
nonautonomous differentiable dynamical systems, and random dynamical systems. 
These vectors identify spatially varying directions of specific asymptotic growth 
rates and obey equivariance principles. In recent years new computational methods 
for approximating Oseledets vectors have been developed, motivated by increasing 
model complexity and greater demands for accuracy. In this numerical study we in- 
troduce two new approaches based on singular value decomposition and exponential 
dichotomies and comparatively review and improve two recent popular approaches 
of Ginelli et al. [17] and Wolfe and Samelson [34J. We compare the performance of 
the four approaches via three case studies with very different dynamics in terms of 
symmetry, spectral separation, and dimension. We also investigate which methods 
perform well with limited data. 

1 Introduction 

The asymptotic behaviour of a linear ODE x(t) = Ax(t), x(t) G M. d is completely deter- 
mined by the spectral properties of the dxd matrix A. Similarly, the long-term behaviour 
of a nonlinear ODE x(t) = f(x(t)) in a small neighbourhood of a fixed point xq, for which 
f(xo) = 0, is completely determined by the spectral properties of the linearisation of / 
at Xq. Well-known extensions of these facts can be constructed when xq is periodic via 
Floquet theory. However, for general time-dependent linear ODEs x(t) = A(t)x(t), the 
eigenvalues of A(t) contain no useful information about the asymptotic behaviour as the 
simple example of P, p. 30] illustrates. On the other hand, if the A(t) are generated 
by a process with well-defined statistics, there is a good spectral theory for the system 
x(t) = A(t)x(t), and this is the content of the celebrated Oseledets Multiplicative Ergodic 
Theorem (MET) ([23], see also Arnold [1] for a thorough treatment), which we state and 
explain shortly. The "well-defined statistics" are often generated by some underlying 
(typically ergodic) dynamical system. 

For clarity of exposition, we will discuss discrete-time dynamics; it is trivial to convert 
a continuous-time system to a discrete-time system by creating eg. time-1 maps flowing 



from time t to time t + 1. Let X denote our base space, the space on which the underlying 
process that controls the time-dependence of the matrices A occurs. As we will place a 
probability measure on X, we formally need a a-algebra X of sets that we can measure^] 
We denote the underlying process on X by T : X O and assume that T is invertible. One 
formally requires that T is measurable^] with respect to X. The "well-defined statistics" 
are captured by a T -invariant probability measure fi on X; that is, \x = \x o T _1 , and we 
say that T preserves /i. Finally, it is common to assume that the underlying process is 
ergodic, which means that any subsets X' G X of X that are invariant (T~ 1 (X') = X', 
implying that trajectories beginning in X' stay in X' forever in forward and backward 
time) have either /i-measure (they are trivial), or /x- measure 1 (up to sets of //-measure 
they are all of X). 

Now we come to the matrices A, which are generated by a measurable matrix-valued 
function A : X — >■ Md(R), where M^(R) is the space of al x d real matrices. We choose 
an initial x G X and begin iterating T to produce an orbit x, Tx, T 2 x, .... Concurrently, 
we multiply ■ ■ ■ A(T 2 x) ■ A(Tx) ■ A(x), and we are interested in the asymptotic behaviour 
of this matrix product. In particular, we are interested in (i) the growth rates 

X(x, v) := lim - log \\A{T n ~ l x) ■ ■ ■ A{Tx) ■ A{x)v\\ 

as v varies in IR d and (ii) the subspaces W(x) C M d on which the various growth rates 
occur. Throughout, ||-|| denotes the standard Euclidean vector norm or the associated 
matrix operator norm \\A\\ = max|| w || = i ||Au||; whether ||-|| is a vector or matrix norm 
will be clear from the context. Surprisingly, the "well-defined statistics" and ergodicity 
ensures that these limits exist, and that there are at most d different values Ai > A2 > 
• • • > > —00 that X(x,v) can take, as v varies over IR d and x varies over //-almost 
all of X (note we allow A^ = —00 to include the case of non-invertible A). We can also 
decompose W 1 pointwise in X as M. d = @f =1 Wi(x), where for all v G Wi(x) \ {0}, one 
has 

lim - log ||yl(T n_1 x) • • • A(x)v\\ = A,. 

The subspaces Wi are equivariant (or covariant) with respect to A over T; that is, they 
satisfy 

W t {Tx) = A(x)Wi(x) 

for 1 < i < i. 

We use the following stronger version of the MET, which guarantees an Oseledets 
splitting even when the matrices A are non-invertible. 

Theorem 1.1 f|14j. Theorem 4.1). Let T be an invertible ergodic measure-preserving 
transformation of the probability space (X,3L,fi). Let A : X — > Md(R) be a measurable 
family of matrices satisfying 

J log + ||A(x)|| d/j(x) < 00. 

1 for example, if X is a topological space, we can set X to be the standard Borel cr-algebra generated 
by open sets on X. 

2 if X is a topological space, and T is continuous, then T is measurable with respect to the standard 
Borcl cr-algebra generated by open sets. 
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Then there exist Xi > X 2 > ■ ■ ■ > Xe > — oo and dimensions mi, . . . , rri£ with m x + • • • + 
me = d, and a measurable family of subspaces Wi{x) C M. d such that for fj,- almost every 
x G X, the following hold. 

1. dim Wi(x) = rrii, 

2. M d = ©f =1 ^(a;) ; 

3. A(x)Wi(x) C WiiTx) with equality if Aj > — oo, 
4- For all v G Wi(x) \ {0}, one has 

lim - log WAiT^x) ■ ■ ■ A(x)v\\ = A*. 

The range of applications of the MET to the analysis of dynamical systems is vast. 
Below, we mention just of few of the settings in which the MET is used. 

Example 1.2. 

1. Differentiable dynamics: One of the first applications of the MET was to differ- 
entiable dynamical systems T : X O on smooth (^-dimensional compact manifolds. 
The matrix function A is the spatial derivative of T, denoted DT. The space M. d 
is associated with the tangent space of X and the equivariance condition becomes 
W t (Tx) = DT(x) ■ Wi(x). If T is uniformly hyperbolic, i:A . >o Wi(x) = W u (x), 
the unstable subspace at x G X and @i.\. <Q Wi(x) = W s (x), the stable subspace 
at x. The spaces Wi(x) provide a refinement of W u (x) and W s (x) into subspaces 
with different growth rates. 

2. Hard disk system: Consider a fixed number iV of hard disks in a region L x x L y 
moving freely between collisions. In each collision a pair of disks change their 
velocities [22]. The region may be finite (hard walls) or periodic (toroidal) in 
either coordinate direction. The quasi-one-dimensional system studied here is a 
two-dimensional system with L y less than twice the particle diameter so that the 
disks remain ordered in the x direction. Here X = ([0, L x ] x [0,^])^ x M? N (with 
the appropriate equivalence classes depending on the choice of hard wall or toroidal 
boundary conditions) is the collection of 4iV-tuples containing all the coordinates 
and momenta of the iV particles. 

The map T : X — > X, x i— > C o J I " r ^ x \x) is the composition of a free-flow map J 7 ^ 1 ') 
and a collision map C. The free-flow map moves the disks in straight lines according 
to their momentum while none of the disks are colliding. The time between collisions 
is the free-flow time t{x) which depends on the initial condition x G X. Collisions 
occur when the boundary of two disks (or one disk and a wall) touch, and the 
collision map exchanges velocities along the direction of collision (since all disks are 
of equal mass). Again, the matrix function A is the spatial derivative of T, so that 
A(x) = DT(x) = D(Co J-^W) (ar). Precise details may be found in [5]. 

3. PDE: The Kuramoto-Sivashinski equation is a model for weakly turbulent fluids 
and flame fronts 
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where v is a damping coefficient. Another familiar example is the complex Ginzburg- 
Landau equation 

r]t = 7} - (1 + i(3)\v\ 2 V + (1 + ia)Vxx 

where T](x,t) is complex and a and (3 are parameters. In both of these cases it is 
possible to approximate solutions of the partial differential equations using Fourier 
spectral methods (see [8 J for details). For instance, in the case of the 1-dimensional 
Kuramoto-Sivashinski PDE we look for solutions of the form 



T)(x,t) = 22 a k(t) e 



ikx/L 



where L is a unitless length parameter, then solve the following system of ODEs 
for the Fourier coefficients a,k(t): 



«fc = [Qk - 1k) a k - «y ^ 



Q'mQ'k—mi 



where qk = k/L. Since the decrease rapidly with k, truncating the above system 
of ODEs is justified. 

In the setting of this review we treat X as the space of Fourier coefficients (ax, ... , a^) 
of the truncated PDE, and consider the transformation T : X — > X defined by 
choosing some r > and letting T ((ai, . . . , a^)) = (ai(r), . . . , Od(r)) where the 
cik(t) are solutions to the system of ODEs with initial conditions 6^(0) = afc. The 
matrix function A is again the spatial derivative of T so that 



da\(r) da\(r) 



A(x) = DT{x) 



8ai 

9Q2(t) 
da i 



9Q2(t) 
da,2 



4. Nonautonomous ODEs and transfer operators: Consider an autonomous 
ODE x(t) = f(x(t)) on X (for example, the Lorenz flow on X = IR 3 ), and its flow 
map £ (r, x) which flows the points x forward r time units. We think of the coordi- 
nates x as a "generalised time" and the ODE x(t) = f(x(t)) is our base system. We 
use this base ODE (the driving system) to construct a nonautonomous ODE or skew 
product ODE as y(t) = F(£(t,x),y(t)). Given an initial time t and a flow time r, 
one may construct finite-rank approximations Px T \t) of the Perron- Frobenius op- 
erator V^ T \x{t)) that track the evolution of densities from base "time" x(t) to 
x(t+r); see [15] for details. The matrices Px T \t) form a cocycle and Oseledets sub- 
space computations enable the extraction of coherent sets in the nonautonomous 
flow (see [IS]). Coherent sets are time-dependent analogues of almost-invariant 
sets for autonomous systems; see [HI [EU H2] . Finite-time constructions for coherent 
sets are described in [16]. In the setting of this review, T : X — > X is defined as 
T(x) = £(t,x), and A{x) = Pt\t). 
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From now on, we denote A(T n 1 x) ■ ■ ■ A(x) as A(x,n). The proof of the classical 
MET [23] proves that the matrix limit 

= lim ({A{x,n))* A(x,n)) 1/2n (1) 

n— toe 

exists for //-almost all The matrix \fr(x) is symmetric, depends measurably on x, 

and its eigenvalues are e Xl ^ > ■ ■ ■ > e Xe ^ x \ The corresponding eigenspaces are denoted 
Ui(x), . . . , Ui(x) and one has 

t i 

Y>[x) , 0u- (,) 0r ; (,). 

Thus Vi(x) captures growth rates from Xi(x) down to \i(x); the "[/"decomposition of 
Vi(x) is orthogonal, while the "W" decomposition (the Oseledets splitting) is equivariant 
(or covariant). 

An alternative notion of stability for non-autonomous systems is the so called Sacker- 
Sell spectrum, cf. [2E]- It is based on exponential dichotomies, cf. [7JI2I] which we briefly 
introduce for linear difference equations of the form 



w n+ \ = A n w n , n G Z, A n G Mrf(R) invertible. (2) 

In the current context we associate the sequence of matrices {A n } n( zz with an invertible 
matrix cocycle over a single orbit, e.g. for some x G X let A n = A(T n x). We restrict 
the introduction of exponential dichotomies to invertible systems only, and note that a 
justification of our algorithm for computing dichotomy projectors strongly depends on 
this assumption. Theory defining exponential dichotomies for non-invert ible matrices is 



contained in eg. [2]. Numerical experiments indicate that Algorithms 3.1 and 3.2 also 
apply in the non-invertible case, however, the corresponding analysis is a topic of future 
research. 

We denote by $ the solution operator of Q, defined as 

A„_i . . . A m , for n> m, 
$(n, m) := ^ I, for n = m, 

An 1 ■ ■ ■ A m-n for n<m. 

Definition 1.3. The linear difference equation ^ has an exponential dichotomy 

with data (K, a s , a u , P*, P%) on J C Z, if there exist two families of projectors P^ and 
P% = I — Pn and constants K, a s , a u > 0, such that the following statements hold: 

P^$(n, m) = $(n, m)P^ Vn, m G J, (3) 



||$(n,m)P^|| < Ke~ a ^ n - m ^ 

Vn > m, n,m G J. 

||$(m,n)P^|| < Ke~ a ^ n ~ m ^ 
Consider the scaled equation 

w n+1 = e' x A n w n , n G Z. (4) 
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Definition 1.4. The Sacker-Sell or dichotomy spectrum is defined as 

°"ed := {A G IR : Q has no exponential dichotomy on Z}. 
The complementary set 1R \ o"ed is called the resolvent set. 

The Sacker-Sell spectrum consists of at most d disjoint, closed intervals, where d 
denotes the dimension of the space, cf. [2H], i-e. there exists an I < d such that 

i 

^ed = [J [AT", A+], where A^ < A7 for i = 1, ...J- 1. 

8=1 

It is well known that the Lyapunov spectrum, when it exists, is a subset of the Sacker- 
Sell spectrum, see [TU] . While the Lyapunov spectrum provides information on bounded 
solutions of $(n, 0), n > 0, the Sacker-Sell spectrum answers this question for $(ra,m), 
n > m. These answers may be different for different initial n because, in contrast to the 
MET setting, there is no a priori stationarity assumption on a base dynamical system 
generating the matrix cocycle. Note that for A Gl \ o"ed h follows from [23J Lemma 2.7] 
that the inhomogeneous equation w n+ i = e~ x A n w n + r n has for every bounded sequence 
r z a unique bounded solution on Z. 

Dichotomy projectors of the scaled equation ^ are constant in resolvent intervals 
Ri : = (Xf, A^L-J, i = 1, ...,£ + 1, where Aq = oo and = -co, see Figure [lj We 
denote these families of projectors by (-P^' s , -P^'")- 

Ri 





CED 




0"ED 


) 









Figure 1: Spectral setup. 
In analogy to the MET we obtain the family of subspaces 

w^ = n(p^ s )nn(p^ u ), nez, i = i,...,e 

that decompose M d for each n e Z 

i=l 

and using the cocycle property ^ it follows for alH = 1, . . . ,1 that 

AnWl = W^ lt neZ. 
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Furthermore, for each w G W^, there exists a constant K = K{w) > such that the 
following equations hold 

||<&(n, ra)iu|| = Ke( Xi +Ti ( n ^ m ))(«- m ) ; f or n > m? w here lim sup rf(n) = 0, 

m)w|| = ife( A ' +Ti ( n - m ))( n_m ) ; f or n <m. where lim sup rj~ (n) = 0. 

n—¥oo 

When working with data over a finite time interval, one has access only to a fi- 
nite sequence A ,A 1 , . . . , Ai-i- 111 this case, one either assumes there is an underlying 
ergodic process generating the sequence A , Ai, . . . , A n _i or one considers exponential 
dichotomies. 

An outline of the paper is as follows. In Sections 2 and 3, we introduce two new 
methods for computing Oseledets vectors. The first method is based on the proof of the 
generalised MET in [TJ] and is particularly simple to implement and fast to execute. 
The second method is an adaptation of an approach to compute dichotomy projectors 
[T5] . In Section 4 we review the approaches by Ginelli et al [T7| and Wolfe and Samelson 
[34] . In Sections 2, 3, and 4 we provide MATLAB code snippets to implement the 
algorithms presented. Section 5 contains numerical comparisons of the performance of 
the four methods on three dynamical systems. The first case study is a dynamical systems 
formed via composition of a sequence of 8 x 8 matrices constructed so that all Oseledets 
vectors are known at time 0; we thus compare the accuracy of the methods exactly in 
this case study. The second case study is an eight-dimensional system generated by two 
hard disks in a quasi-one-dimensional box. The third case study is a nonlinear model 
of time-dependent fluid flow in a cylinder; the matrices are generated by finite-rank 
approximations of the corresponding time-dependent transfer operators. The three case 
studies have been chosen to represent a cross-section of a variety of features of systems 
that either help or hinder the computation of Oseledets vectors, and we draw out the 
advantages and disadvantages of each of the four methods considered. 



2 An SVD-based approach 

The approach outlined in this section is simple to execute and exhibits quick convergence. 
However, as the length of the sample orbit becomes too large this approach fails. 
In [HI proof of Theorem 4.1] it is proven that the limit 

lim A(T- N x,N)Ui(T- N x) 

N— 5>oo 

exists and is equal to the ith Oseledets subspace Wi(x). That is, if one computes [/, in 
the far past and pushes forward to the present, the result is a subspace close to Wi(x). 
Thus, the strategy in [TJ] is to first estimate Z7f in the past and push forward. 

The numerical method of approximating Wj(x), x G X, is implemented in the follow- 
ing steps: 

Algorithm 2.1 (To estimate Wj(x)). 
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1. Choose M, N > and form the matrix 



^(M) ( T -n x ) = (A(T~ N x, M)*A(T- N x, M)) 1/2M (5) 
as an approximation of ([!]) at T~ N x £ X. 



2. Compute U^ M \T n x), the jth orthonormal eigenspace of \I/( M )(T N x) as an ap- 
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proximation of Uj(T x) 



3. Define w\ M ' N \x) = A(T N x,N)Uj(T N x), approximating the Oseledets sub- 



space Wj(x). 

Listing [T] shows part of a MATLAB implementation of Algorithm 2A_ The array 
A = \ [ A(T~ N x) | A(T~ N+1 x) \"' \ A(T M ~ l x)\ contains the d x d matrices which generate 
the cocycle A : X x Z + — > M^R), and the matrix Psi is formed by multiplying the 
matrices contained in A. Step [T] of Algorithm 2.1 is performed prior to the code in Listing 
[TJ Step [2] is performed in lines 1-3 and lines 4-7 perform Step [3j The function returns Wj 
as its estimate to Wj(x). 



Listing 1: Sample MATLAB code of Algorithm 2.1 to approximate WAx 



svd (Psi) ; 



sort(diag(s) , 'descend ' ) ; 



1 [ ,3, 

2 [~,P] 

3 Wj = v ( : , p( j) ) I norm ( v ( : , p( j ) ) ) ; 

4 for h = 1 : N 

5 Wj = A(: , (h—1)* dim+l:h* dim)* Wj ; 

6 f/7 = f/j'/ norm (Wj); 

7 END 



The values of M and iV can be chosen with relative freedom and in our examples that 
follow we have chosen M = 2N to compute over a time window centred on x, from T~ N x 
to T N x. Unfortunately, we cannot choose M and N arbitrarily large and expect accurate 
results. If A(T~ N x, M) is constructed via the product 

A(T~ N x, M) = A(T M ~ N ~ 1 x) ■ ■ ■ A(T- N+l x)A(T~ N x) 

then with larger M the numerical inaccuracies of matrix multiplication compound and 
this product becomes more singular and thus a poorer approximation of A(T~ N x,M). 
Because of this, ^^ M \T~ N x) cannot be expected to accurately approximate ^(T~ N x) for 
large M. However, even if we suppose ty( M ^(T~ N x) accurately approximates ^(T~ N x), 
the small, but non-zero, difference in Uj M \T~ N x) and Uj(T~ N x) grows roughly as 
O (e Ar ( Al-Aj ' ) ) during the push-forward in step [3] above. For these reasons M and iV 
must be chosen carefully. 



2.1 Improving the basic SVD-based approach 

We present a simple improvement that can overcome one of the sources of numerical 
instability, namely the push-forward process in step [3] above. 
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Figure 2: Schematic of the re-orthogonalisation described in Section 2A_ The black line 
represents the orbit centred at x G X and the points T~ Nk x are those points at which 

To do this we use the (blue) 
x) 1 - and perform the (red) push- 



-N k 



we ensure orthogonality with the subspaces Vj(T 
approximations ^^ M \T~ Nk x) to approximate Vj(T~ 



x ) 



forward and orthoganlisation steps starting with [/■ (T Nl x) and ending with Wj 



r(Af,0). 



X 



(see Algorithm 2.2) 



Recall that the subspace Vj(x) = Uj(x) © • • • © Ut(x) = {U\{x) © • • • © Uj-±(x)) ± is 
A-invariant and that for v G Vj(x)\Vj+i(x) (with Ve + i(x) = {0}) we have 

Xj(x) — lim — log \\A(x, n)v II . 

71— >OD n 

The subspace Vj(x) contains Wj(x), Wj+i(x), . . . , We(x), and so the Oseledets sub- 
space Wj(x) is necessarily perpendicular to all Ui(x), . . . , Uj-i(x). To solve the numerical 
instability of step [3] we enforce this condition periodically. 

The amended algorithm is implemented as follows: 

Algorithm 2.2 (To estimate Wj(x)). 

1. Choose M, Ni > N2 > ■ ■ ■ > N n = and form the matrices 

¥ M \T- Nk x) = (A(T- Nk x,M)*A(T~ Nk x,M)) 1/2M , k = l,...,n. 

2. Compute all the orthonormal eigenspaces U^ M \T~ Nk x), i = — 1 of (5p 
(replacing N with in (5)) and the eigenspace Uj (T~ Nl x). 

3. Let projy : M d — > M. d be the orthogonal projection onto the subspace V so 

that ker (proj v ) = V 1 and vj M \x) = (u[ M \x) © • • • © U^(x)\ . Define 

Wj' Nl \T~ Nl x) = Uj (T~ Nl x), and then define iteratively by pushing forward 
and taking orthogonal projections: 

w (M, Nk+l) {T -N k+lx) = W0 ) vr{T - Nk+lx) (A (T- N *x,N k+1 - N k ) W}"™ (T~ Nk x) 

r(M,N n ) / \ _ rrAMfi) 



4. Wj ' " (x) = Wj ' (x) is our approximation of Wj(x). 



j 
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Listing [2] shows an example implementation of Algorithm 2.2 in MATLAB. Lines 1 



18 are responsible for performing Steps [T] and [2j whilst the push forward procedure of 
Stepplis performed in lines 20-30. Again, the matrix cocycle is stored in A = \A(T~ N x) | 
A(T x) I • • • I A(T M ~ 1 x)\ and the function returns Wj as its approximation of Wj(x). 
The variable Nk is a one- dimensional array containing the elements of {Nk} and is counted 
by k. 



Listing 2: Sample MATLAB code of Algorithm 2.2 to approximate Wj(x 



1 fc=0; 

2 for n = iVfc, 

3 Psi = eye ( dim) ; 

4 for i—0:M— 1, 

5 Psi = A ( : , ( n + i — l)*di/?? + l:(n + i)* dim)* Psi ; 

6 Psi = Psi/NORMEST(Psi ) ; 

7 END 

8 [~ , s , u] = SVD (Psi ) ; 

9 [~ , jj] = sort ( di AG ( s) ,' descend ') ; 

10 if n==l, 

11 Wj = u(: ,p( j))/norm(u(: 

12 ELSE 

13 FOR i = 1 : j7" — 1 , 

14 fc = fc+1; 

15 U(:,i,k) = u(: ,p( i) ) / norm ( u ( : i) ) ) ; 

16 END 

17 END 

18 END 

19 fc=0; 

20 for n = l:N, 

21 Wj = A(: , ( n—l)* dim+l:n* dim)* Wj ; 

22 f/7 = {/ji'/nqrm( f/j/') ; 

23 if ANY (Nk = n+1), 

24 fc = fc+1; 

25 FOR i = 1 : J — 1, 

26 Wj = Wj - DOT ( Wj , CT( : , i , fc) ) * U( : , i , fc) ; 

27 f/j = Wj/ NORM ({/?); 

28 END 

29 END 

30 END 



Remark 1. Unfortunately, some numerical issues with this approach remain. They stem 
primarily from the long multiplication involved in building the variable Psi of Listings [T] 



and 2 This results in Psi becoming too singular and hence uj- (T n x) (j 7^ 1) poorly 
approximates Uj(T~ n x). As can be seen in Section [5j Algorithm 2.2 works superbly for 
~ n x) is well approximated for large n. However when estimating Wj(x), 
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j > 2, a good estimate of U^f}{x) is required for an accurate projection proj y (M) and for 
large n such an estimate becomes unreliable. 

3 A dichotomy projector approach 

We derive an approach for the computation of a vector w{ G W-{ = 1Z(P^ S ) D lZ(P-{ +1,u ). 
For this task, we first need a guess of A nght G Ri and of A left G -R«+i in two neighbouring 
resolvent intervals that lie close to the common spectral interval, see Figure [3] 



^lcft fright 



t t 



Ai 
+ 



Rs R2 
Figure 3: Choice of A left / right in case i — 2. 



Numerical experiments indicate that we get the best results by choosing A nght and 
A left close to (but outside) the second Sacker-Sell interval. This conclusion is supported 



by theoretical estimates on the approximation error for Algorithm 3.2, discussed at the 
end of Section 3. 

The following observation from [HI [19] allows the computation of dichotomy projec- 
tors by solving 

w l n+1 = A n w % n + 8 n ^ m -iei, n G Z, e; i-th unit vector. (6) 

With Green's function, cf. [24J, the unique bounded solution w l z of ^ has the explicit 
form 

< = G{n, m)e h n G Z, where G{n, m) = I m)P " 1 ' " " m ' (71 



and consequently 

p s = I w l ... w a 



rn 



Numerically, we approximate the unique bounded solution on Z by the least squares 
solution of ^ on a sufficiently long interval. For an error analysis of this approximation 
process, we refer to [HI Theorem 4]. 

The algorithms that we propose in this section compute a vector w G Wq in analogy 
to Wj(x) in the previous sections. For simplicity, we restrict the representation to the 
case j = 2 and assume that and W% are one-dimensional subspaces. 
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In the absence of information about the dichotomy intervals, one may proceed as 
follows. Given a finite sequence of matrices, one can estimate a point in the spectral 
interval [\~X£],q = 1,2,3 by computing the (logarithmic) growth rates of one-, two-, 
and three-dimensional subspaces using direct multiplication; these growth rates should 
approximate Ai, Ai + A 2 , and A1 + A2 + A3, respectively. By taking differences to obtain esti- 
mates X q , q = 1,2,3, (the caret indicating estimated quantities) one should obtain values 
in the interior of [A~ , A+],g = 1,2,3. We then estimate A lcft = A 2 — (A 2 — A 3 )/10 ;$ A 2 

and A ri § ht = A 2 - (A 2 - A^/IO > A+. 

In the first step of our first algorithm, we compute a basis of the two-dimensional 
space 1Z(Pq ,u ). Then, in the second step, we search for the direction w in this subspace 
that additionally lies in 7Z(P ' S ) and assure in this way that w G n(P 2 ' s )nn(P^ u ) = W 2 . 

Algorithm 3.1 (A Dichotomy Projector approach to estimate W 2 (x) by computing Wq). 



1. Suppose N eN and consider n G [-N, N] n Z. Let A n = A(T n x). 
Solve the least squares problem 



n+l 



A n w l n + 8 n ^r\ n = -N, . . . , N - 1, i = 1, 2 
such that IKwljy, . . . ,w]v)ll2 is minimised, 



(8) 



where the r l are chosen at random and ||-|| 2 is the ^ 2 -norm. Define p 1 := A^iW^, 
i = 1, 2. 



2. Solve for W[o,jv] and n the least squares problem 

w n+1 = e ~ Arisht A n w n , n = 0, . . . , N - 1, 

W + Kp 1 + p 2 = 0, 

such that \\(wq, . . . , wn, k)\\ 2 is minimised. 
Then w is our approximation of w 2 (x) G W 2 {x). 



(9) 
(10) 



The unique bounded solutions on Z of these two steps satisfy p 1 ,p 2 G TZ(Pq U ) and 
these vectors are generically linear independent. Furthermore w G 1Z(Pq' s ) due to (|9| 
and w G TZ{P^ U ) due to @. Thus w G TZ{P^' U ) n n(P 2 ' s ) = W 2 . 

Note that (|8l) has the form 



where 



Bw = r, with B G M 2dNA2N+1) (R), r G R 2dN , 

I \ ( w^ N \ 



B 



f-e 



-A lett A 



—e 



w 



A N I) 



\W N -1 J 



and the nth entry of r is the vector 5 n -ir l G M. d for i — 1, 2. 

The least squares solution can be obtained, using the Moore-Penrose inverse: 



w = B + r, where B + = B T (BB 



T\-l 
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Listing 3: Sample MATLAB code for Algorithm 



3.1 



1 % step 1 

2 B = zeros (2*N* dim, 2* ( #+1)* dim) ; 

3 for i = 1:2* N 

4 B( dim*(i — 1) + 1: dim* i , dim*(i — l) + l: dim* i) 

5 = — exp(— Lambda_ left)*A(: , dim*( i — 1)+1: dim* i ) ; 

6 B(dim*(i — 1) + 1: dim* i , dim*( i) + 1: dim*( i + 1)) 

7 = EYE ( dim) ; 

8 END 

9 R = zeros (2* dim*N, 2 ) ; 

10 R( dim*(N— 1)+1: dim*N, : ) = rand ( dim , 2 ) ; 

11 y = (B*B')\iJ; 

12 u = B' * y; 

13 p.Z = A ( : , dim*(N— 1)+1: dim*N) *u( dim*(N— 1) + 1: dim*N, 1 ) ; 

14 p2 = j4 ( : , dim*(N— 1) + 1: dim*N)*u( dim*(N—l) + l: dim*N, 2 ) ; 

15 pi = t'J/NORM ( t/J ) ; 1/2 = f2/N0RM ( i/2) ; 

16 % siep £ 

17 £ = zeros ( dim*(^+l) , dim*(#+l) + l); 

18 for i = 0:#-l 

19 B( dim* i + 1: dim* (i + 1) , dim* i + 1: dim* (i + 1)) 

20 = — exp(— Lambda_ right)* A (: , dim*( i+N) + l: dim*( i+N+1)) ; 

21 B( dim* i + 1: dim* ( i + 1) , dim*( i + l) + l: dim* ( i + 2)) 

22 = EYE ( dim) ; 

23 END 

24 B( dim*N+l: dim*(N+l) ,1: dim) = EYE(dim); 

25 B(dim*N+l:dim*(N+l) ,dim*(N+l)+l) = pi; 

26 R = ZEROS ( dim* (N+l ), 1) ; 

27 R( dim* N+l: dim* (N+l) ,1) = -p2; 

28 y = (B*B')\R; 

29 u = B'*y; 

30 iu2 = u( di/w*#+l: dim*(#+l))/NORM («( dim*N+l: dim*(N+l)) ) ; 



cf. [5D]. Numerically, we find iD by solving the linear system BB T y = r; then w = B T y. 
Note that in the unlikely case where p 1 G Wq, Algorithm 3.1 fails. An alternative 



approach for computing vectors in Wq that avoids this problem is introduced in Algorithm 
3.2| The main idea of this algorithm is to take a random vector r, project it first to 



1Z(Pq ,u ) and then eliminate components in the wrong subspaces, by projecting with Pq' s . 

Algorithm 3.2 (An alternate Dichotomy Projector approach). 

1. Again, suppose N e N and consider n E [— N, N] fl Z and let A n = A(T n x) as 
above. Solve the least squares problem 

w n+1 = e~ AlcCt A n w n + 5 n _ir, n = -N, . . . , N - 1, (11) 
such that ||(w_iv, . . . ,w)jv)||2 is minimised, 
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where r is chosen at random, and define r' = 
2. Solve the least squares problem 



w' 



wl 



n+1 = e - Arisht A n w' n + S^-tr', n = -N, . . . , N - 1, (12) 
such that IKw^jv, . . . , i^jv) 1 1 2 is minimised. 

Then w' is our approximation of w 2 (x) £ W 2 {x). 

The solution wq on Z of these two steps satisfies wq = P^P^r £ 1Z(Pq ,s )P\1Z(Pq' u ) = 



Listing 4: Sample MATLAB code for the second step of Algorithm 3.2 



IB — zeros ( 2 * TV* dim, 2* ( #+1)* dim) ; 

2 for i = 1 : 2 * iV 

3 B( dim*(i — l) + l: dim*i , dim*(i — l) + l: dim* i) 

4 = — exp(— Lambda_right)* A(: , di/??*( i — 1)+1: dim* i) ; 

5 B( dim* ( i — 1)+1: dim* i , dim* i+1: dim* ( i + 1)) = EYE(dim); 

6 END 

7 R = zeros (2* dim*N } 1 ) ; 

8 R(dim*(N-l) + l:dim*N,:) = pi; 

9 y = (B*B')\R; 

10 14 = 

11 w2 = u( dim*#+l: dim* (#+1))/norm ( u( dim* N+1: dim* (N+ 1 ) ) ) 



3.1 Error estimate 



We give an error estimate for the solution of Algorithm [3^ for a finite choice of N. Details 
on deriving this estimate are postponed to a forthcoming publication. 

For A and A right close to the boundary of the second Sacker-Sell spectral interval, 
we denote the dichotomy rates of 

A left . A right , v 

w n+ i = e A n w n , w n+1 = e A n w n , n £ Z (13) 



by (a ,s ,a ,u ) and (a r,s , a r ' u ) , respectively. Let w be the solution of Algorithm 3.2 on Z 
and let Wq be its approximation for a finite choice of N. Careful estimates show that the 
approximation error in the "wrong subspace" 1Z(Q), with Q :— I — Pq' s Pq' u is given as 



\\Q(w -w )\\ <CN(e-^ SN + e- ar ' UN ), (14) 

where the constant C > does not depend on N. 

The exponential dichotomy rates a e ' s and a r,u of the difference equations ( 13 ) depend 
on the choice of A left and A right in the following way: for A left in the resolvent set -R3 = 
[Aj!~, A 2 ] the difference equation 

A left . 

w n+1 = e A n w n , n £ Z 
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has an exponential dichotomy with stable dichotomy rate a e ' s for all a e,s with 

< a e ' s < A lcft - A+ 
Similarly, for A nght in the resolvent set R2 = [Aj, A^] the difference equation 

-A right A ^ H7 

w n+ i = e A n w n , n e Z 

has an exponential dichotomy with unstable dichotomy rate a r,u for all a r ' u with 

< a r ' u < Xi - A right . 

Note that both of the above inequalities are strict. 

Inspecting equation (14), we get the best (smallest) maximal error if we choose A 
i?3 and A nght e R 2 so as to maximise a e,s and a r ' u . Consequently, we get the best 
numerical approximations, if A left and A nght are chosen close to, but not equal to, the 
boundary of the common spectral interval [A^", Aj]. 



left 



4 The Ginelli and Wolfe schemes 



4.1 The Ginelli scheme 

The Ginelli Scheme was first presented by Ginelli et al. in [17J as a method for accu- 
rately computing the covariant Lyapunov vectors of an orbit of an invertible differentiable 
dynamical system where the A(x) = DT(x) are the Jacobian matrices of the flow or map. 

Estimates of the Wj(x) are found by constructing equivariant subspaces Sj(x) = 
Wi(x) © ■ ■ ■ © Wj(x) and filtering the invariant directions contained therein using a power 
method on the inverse system restricted to the subspaces Sj(x). 

To construct the subspaces Sj(x) we utilise the notion of the stationary Lyapunov 



basis [UJ. Choose j orthonormal vectors si(T n x),S2(T 
that Si(T~ n x) ^ Vj + i(T~ n x) for 1 < i < j and construct 



x) 



Sj(T n x), n > 1, such 



~(n) 
S) (X 



(x) = A(T- n x,n) Sl (T- n x), i 



■ J- 



Using the Gram-Schmidt procedure, construct 
{si(x), . . . , Sj\x)} from {s[ n ^(x), . . . , s™ (x)}, that is, 



the orthonormal basis 



4 n) (x) 



s[ n \x) 



4 B) (x)-(4 B) (x). 8 W(x)) sf\x)) 



(15) 
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Then as n — > oo the basis {s^ (x), . . . , Sj(x)} converges to a set of orthonormal vectors 

{si (x), . . . , Sj \x)} which span the j fastest expanding directions of the cocycle A [TT], 
that is, if the multiplicities mi — • ■ ■ — rrij — 1 



Sjix) := span {s^°\x) } . . . } sf°\x)^ = Wi(x) © • • • © t^-(x) 

= \/ 1 (x)\V 7+1 (x). (16) 



If the Oseledets subspaces are not all one-dimensional, that is the Lyapunov spectrum 
is degenerate, then we choose Sj(x) only for those j which are the sum of the first k 
multiplicities, i.e., j = mi + ■ ■ ■ + m^. Then 



In the interest of readability we assume the Oseledets subspaces are one-dimensional but 
note that the approach may be extended to the multi-dimensional case. 
Note that the Sj(x) are equivariant by construction: 



provided j < mi -\ + mi-i if \i = — oo. 

We describe the Ginelli approach to finding Wv(x). Suppose dimWi(x) = 
dimW / 2(x) = 1 and Ai > A2 > —00 and that the basis {s[ (x),s^ (x)} is known 
at x E X. Note first that span |si°°^(x)| = Wi(x). Let c(x) E M 2 denote the coefficients 

of w 2 (x) E W 2 (x) in the basis {s ( f°\x), s 2 °°\x)} (recall that the orthogonal projection of 
w 2 (x) onto s[°°\x) is zero for i — 3, 4, . . . , d) 



Lemma 4.1. Let Q(x) denote the d x 2 matrix whose ith column is s 4 (x). Then for 
each n > there exists an upper triangular, 2x2 matrix R(x, n) satisfying 



S ] (x) = W l (x)@---®W k (x) 
= Vi{x)\V k+ i{x). 



A(x,n)Sj(x) = Sj{T n x) 



then 



w 2 (x) = ci(x)si'(x) + c 2 (x)s 2 '(x). 



A{x,n)Q(x) = Q(T n x)R(x,n). 



(17) 



Proof. Note that 



A(x,n)Q(x) = A(x,n) sf^x) s^\x) 





Q(T n x)R(x,n) 
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where 



Q(T n x) 



(T n x) 



(T n x) 



and 



R(x, n) 



A(x, n)^ (x)\\ (s^ (T n x) , A(x, n)4°° ) (x) 



using the equivariance of = span{s^ (a;)} and S2 (x) = spaD.{sf° > (x), s^'fa)}. □ 



(oo). 



,(o°)/ 



Thus, the QR-decomposition of Lemma 4.1| is equivalent to the Gram-Schmidt or- 
thonormalisation that defines the stationary Lyapunov bases. The columns of Q(T n x) 
form the stationary Lyapunov basis at T n x. 

We have chosen the above notation R(x, n) specifically since, defined in this way, R 
forms a cocycle which is the restriction of A to the invariant subspaces Sj. 

Lemma 4.2. The matrix R(x,n) defined above forms a cocycle over T. 

Proof. Let n, m > then 

A(x, n + m)Q(x) = Q(T n+m x)R{x, n + m) (19) 



by Lemma 4.1 Since A(x, n + m) = A(T n x, m)A(x } n), 

A(x, n + m)Q(x) = A(T n x, m)A(x, n)Q(x) 

= A{T n x,m)Q{T n x)R{x,n) 
= Q(T n+m x)R{T n x,m)R{x,n). 



(20) 



Equating (19) and (20) gives 

R(T n x, m)R(x, n) = R(x, n + m) 
as Q(T n+m x) is left-invertible. 



□ 



Since c(x) is the vector of coefficients of the second Oseledets vector of the cocycle A, 
it is the second Oseledets vector of the cocycle R. To see this, recall w 2 (x) = Q(x)c(x) G 
W 2 (x) so that 



A 5 



which, due to (17), becomes 

A 2 : 



lim —log\\A(x,n)Q(x)c(x) 

n— >oo fl 



lim -log\\Q(T n x)R(x,n)c(x) 

n— >oo fl 

lim — log || R(x, n)c(x) \\ 

n— >oc fl 



since the columns of Q(T n x) are orthonormal. 

We may approximate c(x) numerically using a simple power method on the inverse 
cocycle R^ 1 (which exists since Ai > A2 > —00). 

The Ginelli method can be summarised by the following steps: 
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Algorithm 4.3 (Ginelli method of approximating W2{x) G Wzfa)). 

1. Choose x G X and M > and form {s^(x), s^^)} by first randomly se- 
lecting two orthonormal vectors {si(T~ M x), S2{T~ M x)} then performing the push- 



forward/Gram-Schmidt procedure given by (15). That is, define 



sf\ x ) = A(T- M x,M) Si (T- M x), i = 1,2, 



followed by setting 



4 M) (o;) = M (s?\x) - (4 M) (x) ■ s[ M \x)) s[ M \x) 

where M : v h-> t>/ ||t>||. The vectors {s^ (V), ^^(x)} form an approximation to 
the stationary Lyapunov basis {si\x), (x)}. 



2. Choose N > and using the approximate basis {s[ M \x), s^\x)} in (18), form an 
approximation to R(x,N), denoted by R^ M \x,N). 

3. Choose d G M 2 either at random, or by some guess at the second Oseledets vector 
of R at T N (x) G X, in this review we found d = (0, 1) to work well. Use the inverse 
iteration method to approximate c(x), that is, define our approximation to c(x) as 

c^ N \x) = RW\x,N)- 1 d 

= R( M \T N x,-N)c' 

4. Then 



wi M ' N) (x) = s[ M) (x) s{ M) (x) c^ N Hx) 



is our approximation to w 2 (x) G W 2 (x). 

As before, there is some freedom of choice of both M and N as well as of the initial 
orthonormal basis {si(T~ M x), S2(T~ M x)}, used to approximate S^fx), and of the 2-tuple 
d. The larger M and N are chosen, the more accurate w*^ ,N \x) will be, provided 
s 2 {T- M x) i Wi(T- M x)UV z (T~ M x) and d £ £i(T M x) where E x is the Oseledets subspace 
of R with Lyapunov exponent Ai. 

Listing [5] shows an example implementation of Algorithm 4.3 in MATLAB which 



approximates Wj(x) G Wj(x) using M = N and si(T M x), . . . , Sj(T M x) are chosen 
at random and d = (0, ...,0,1). Lines 2 through 6 construct the approximation of 

j—l entries 

the stationary Lyapunov basis, ^s[ M \x) , . . . , s^ M \x)^ which is stored as columns of 

the matrix Q(x) represented as the variable QO, while lines 7 through 11 construct the 
cocycle R stored in AllR as [R(x) \ R(Tx) | • • • | R(T x)\. Lines 12 through 17 perform 
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Listing 5: Sample MATLAB code for Algorithm 4.3 



1 [ Q , ~ ] = qr(rand ( dim , j) , ) ; c = [ zeros (1 , j— 1) 1 

2 for i — l:N, 

3 QNew = A(: , ( i — 1)* dim+1: i* dim)* Q; 

4 [Q,~] = qr( QNew, 0); 

5 END 

6 QO = Q; 

7 for i = N+l:2*N+l, 

8 QNew = A(: ,(i — l)*dim +1: i* dim)* Q; 

9 [Q,R] = QR(QNew,0) ; 

10 AllR = horzcat (R, AllR) ; 

11 END 

12 numOfR = size( AllR ,2) / j ; 

13 for i = l:numOfR 

14 # = AllR(: , ( i — 1)* j'+l: i* j) ; 

15 cilfew = ii\c; 

16 c = ciVeiy/NORM ( cNew) ; 

17 END 

18 w = #0*c; 



a simple power method on R(x, N)' 1 to find the coefficient vector c, which represents the 
approximation of Wj(x) in the basis js^ (q;), . . . , s^ M \x) |. Thus, the approximation is 
given by Q(x)c. Although Algorithm 4.3 is specific to the case where j = 2, Listing [5] is 



applicable to any j for which R(x, N) 1 exists. 

It can be shown that, in this case where the top Lyapunov exponent has multiplicity 
1, E 1 = span{(l,0) T }. 

Lemma 4.4. // the first Lyapunov exponent has multiplicity 1, the dominant Oseledets 
subspace of the cocycle R is E\ = span{(l, 0) T }. 

Proof. Recall that s < f°\x) G Wi(x) since span js^(:r)j = Si(x) = Vi(x)\V2(x) (from 
and for all s e ^i(z)\F 2 (^) 

Ai = lim — log \\A(x, n)s\\ . 

n->oo n 

We may write s = Q(x)(a, 0) T for some o6K then 

Ai = lim — log II A(x, n)Q(x)(a, 0) 



T\ 



n— >oo Ti 

lim - \og\\Q{T n x)R{x,n){a,0) 



n— too fl 

and since ||Q(T n x)s|| = (the columns of Q are ortho normal) 

1 



Ai = lim — log \\R(x,n)(a, 0) 

n— >oo fl 



Tl 



Since S\ is A-invariant, span{(l,0) } is .R-invariant and the proof is complete. □ 
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4.1.1 Limited Data Scenario 



In the case where convergence isn't satisfactory because the amount of cocycle data 
available is too small (for any M and N to be small), the approximations from Algorithm 



4.3 can be improved by using better guesses at si(T M x), s 2 (T M x) and d . 

Note that s[ (x) and s^°\x) are two orthonormal vectors optimised for maximal 
growth over the time interval [—00, 0] D Z. In the situation where the values of M and 
iV are limited, one can choose those two vectors that are optimised for growth over the 
shorter time interval [— M, 0] n Z. In [31] this is achieved by computing the left singular 
vectors of A(T~ M x, M). This approach works well for small M but can become inaccurate 
for very large M for the reasons in Remark [1} 



In practice, we have observed that a combination of Step 1 in Algorithm 4.3 and the 
above provides the most robust method of accurately approximating s[ (x) and s% (x). 



As an alternative to Algorithm 4^3 the following may be used: For M > M', compute 
vectors optimised for growth from — M to — M + M', then push-forward these vectors 
from —M + M' to 0. 



Algorithm 4.5 ( Improved Algorithm 4.3). 



1. Choose x € X and M > M' > 0. Compute the two left singular vectors 
of A(T~ M x, M') corresponding to the two largest singular values and call them 

~ Sl {T- M+M 'x) and s 2 (T~ M+M 'x). Now define {s[ M>M '\x), s { 2 M ' M '\x)} as an ap- 
proximation to |s^°°^(a;), S2°°^(a;)| by the Gram-Schmidt orthonormalisation of 
A(T~ M + M 'x, M - M') Sl (T- M+M 'x) and A(T~ M+M 'x, M - M')s 2 (T- M+M ' x) as in 



(15). 



Steps 2-4 as in Algorithm 4.3 



In practice, one should choose M' large enough so that enough data is sampled, but 
not so large that A(T~ M x, M') is too singular. 

4.2 The Wolfe scheme 

The approach followed by Wolfe et al. [34j directly computes the subspace splitting 
as the intersection of two sets of invariant subspaces. The description of the numerical 
construction of the subspaces Sj(x) featured below differs slightly from [31], however, 
the essential features of the approach are retained. In fact, the constructions featured 
here improve upon those in [34J in terms of accuracy versus amount of cocycle data used 



- in the notation of Algorithm |4.6| below, if M\ is made larger, w^ 2 Ml ' Ml ' M2 \x) is more 
accurate, which is not the case in [34J for the same reasons discussed in Remark [Tj 

Recall the eigenspace decomposition Uj(x) of the limiting matrix ^>{x) presented in 
Section |2] and define Vj(x) = Uj(x) © ■ ■ ■ © U t (x). Recall that Vj(x) D Wj(x),W j+1 (x), 
. . . , Wt(x). Also recall from the previous section that Sj(x) D W\(x), W 2 (x), . . . , Wj(x). 
Thus 

Wj(x) = Vj(x) n Sj(x). 
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Again, in the interest of readability we assume the Oseledets subspaces Wj(x), and 
the eigenspaces Uj(x), are one-dimensional. As in the case of the previous section, the 
ideas here may be extended to the case in which the Oseledets subspaces are not one- 
dimensional. Let Uj(x) be the singular vector spanning Uj(x) and let Sj(x) = Sj \x) be 
the jth element of the stationary Lyapunov basis as in the previous section. Then note 

3 

i=l 
d 

W j{ X ) = ^2 ( w j( x )> u i( x )) Ui(x). 
i=j 

Taking inner products with u k {x) and s k {x) respectively gives 

3 

(wj(x),u k (x)) = 22 ( w j( x )' s i( x )) (si(x),u k (x)) for k>j, 
i=i 

d 

■ Sk(x)) = "22 ( w j( x )i u i( x )) ( u i( x )i s k(x)) for k < j. 



(21) 
(22) 



Substituting (21) into (22) and rearranging gives 

3/d \ 
(wj(x), s k (x)) = ^2 I ^2 ( s k(x),u h (x)) (uh(x), Si(x)) J (Si(x),Wj(x)) . (23) 

t=l \h=j J 

Note that Yl=i ( s k(x),u h (x)) (u h (x), Si(x)) = 8 H so 

d j-i 
"22 (s k (x),u h (x)) (u h (x), Si(x)) = 5 ki - ^ (s k (x),u h (x)) (u h (x), Si(x)) . 

h=j h=l 



Then (23) becomes 



j j'-i 



(Wj(x),s k (x)) = ^25 ki (si(x),Wj(x)} - ^2^2 ( s k(x),u h (x)) (u h (x),Si(x)) (Si(x),Wj(x)) 



i=l 



i=l h=l 



j J-l 



(s k (x),Wj(x)) - y]y] (s h (x),u h (x)) (u h (x),Si(x)) (Si(x),Wj(x)) 



i=l h=l 



j i-1 



= ^2 ^2 ( s k(x),u h (x)) (u h (x), Si(x)) (Si(x),Wj(x)) . 



(24) 



i=l h=l 



Equation (24) may be considered as a j x j homogeneous linear equation by defining a 
matrix entry-wise as 

j'-i 

D k i = "22 ( s k{x),u h (x)) (u h (x),Si(x)) 

h=l 
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and solving 



Dy = 0, 

where y^ = (si(x),Wj(x)). The entries of y are then the coefficients of Wj(x) with respect 
to the basis sx(x) . . . , Sj(x). 

The Wolfe approach may be implemented as follows: 

Algorithm 4.6 (Improved Wolfe approach to approximating Wj(x) G Wj(x)). 

1. Choose x G X and Mi > M[ > and construct {s[ Ml,Ml \x), . . . , s^ Ml ' Ml \x)} as 
an approximation of the stationary Lyapunov basis vectors {s^ (a;), . . . , s^f°\x)} 



using the methods outlined in Step 1 of Algorithm |4.5[ that is, compute the left 
singular vectors of A(T~ Ml x, M[) corresponding to the j — 1 largest singular values 

and call them Si(T~ Ml+M ^x), i = 1, ... ,j — 1. Then form the s[ Ml ' Ml \x) by the 
Gram-Schmidt orthornormalisation of A(T~ Ml+M tx, Mi-M[)§i(T- Ml+M 'ix) for i = 
1 ./-• 1- 

2. Choose M2 > and construct the one-dimensional eigenspaces u[ M ' 2 \x),..., 
(x) as approximations to the eigenspaces U\(x), . . . ,Uj-i(x) as in Step 1 of 



Algorithm 2.1 that is, construct 

¥ Ah \x) = (A(x 1 M 2 )*A(x,M 2 )) 1/2M2 , 

and let U\ M% \x) be the zth orthonormal eigenspace of ^^ M2 \x). Define u^ l2 \x) G 
U^\x),i = l,...,j-1. 

3. Form the matrix D as above: 



h=l 



4. Solve the homogeneous linear equation Dy = 0. Then w ( Ml ' Ml > M ^ = J2i=iVi s 
forms our approximation of Wj(x) G Wj(x). 



x 



This approach suffers from the same numerical stability issue of Algorithm 2.2 of Sec 




Namely, the vector spaces Lf[ M \x), . . . , U^(x) may only poorly approximate 



, Uj-\(x) for M too large (see the final paragraph of Section 2.1). 
A recent paper [20J provides alternative descriptions of both the Ginelli et al. and 
Wolfe and Samelson methods, well-suited to those familiar with the QR-decomposition 
based numerical method for estimating Lyapunov exponents due to Benettin et al. [31 H] 
and Shimada and Nagashima [29] . The discussion in [20] is restricted to invertible cocycles 
generated by the Jacobian matrices of a dynamical system. Although this assumption 
allows stable numerical methods to be constructed, i.e., better convergence obtained for 
larger data sets, it means some important examples in which the matrix cocycle is non- 



invertible are overlooked, for example the case study of Section 5.3 While the memory 



footprint of the implementations discussed in J2D] is estimated, there is no discussion 
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of convergence rates or accuracy with respect to the amount of cocycle data available. 
Finally, while the examples featured in [20] explain the methods presented in the context 
of differentiable dynamical systems the case studies of Section [5] in the present paper 
focus on comparative performance of the methods presented, via a broad range of possible 
applications. 



5 Numerical comparisons of the four approaches 

We present three detailed case studies, comparing the four approaches for calculating 
Oseledets subspaces. The first case study is a nontrivial model for which we know the 
Oseledets subspaces exactly and can therefore precisely measure the accuracy of the 
methods. The second case study produces a relatively low-dimensional matrix cocycle, 
while the third case study generates a very high- dimensional matrix cocycle; in these case 
studies we use two fundamental properties of Oseledets subspaces to assess the accuracy 
of the four approaches. 



5.1 Case Study 1: An exact model 

In general the Oseledets subspaces cannot be found analytically which makes the task of 
determining the efficacy of the above approaches difficult. However, the exact model de- 
scribed below allows us to compare the numerical approximations with the exact solution 
by building a cocycle in which the subspaces are known a priori. 

We generate a system with simple Lyapunov spectrum Ai > A2 > • ■ • > A^ > —00. 
We form a diagonal matrix 



R 



0\ 







3 Ad 



and generate the cocycle by the sequence of matrices {A n } where 
A n = S n RS n _ 1 

for n ^ -l,n G \-N,N] n Z, 



S „ 



I + eZ, 

/o 



z 2 



V 



\ 



z d J 



for n = — 1. 



(25) 



The entries of Z and the numbers z 2 , ■ ■ ■ , Zd are uniformly randomly generated from the 
interval [0, 1]. By construction, the columns of SV1-1 span the Oseledets subspaces at time 
n E [-N, N] n Z. 

We compare the exact result at time n = with the approximations computed by 
the various algorithms for d — 8, {Ai, . . . , Ag} = {log 8, log 7, log 6, . . . , log 1} and e = 0.1 



for varying amounts of cocycle data {A(T 



N 



X 



A{T N x)}. The exact model has a 
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well separated spectrum, is generated using invertible matrices, and is of relatively low 
dimension. 

For this model we use the following choice of parameters to execute the algorithms: 



Algorithm |2j2| M = N and {N k } = {1, 6, . . . , 5k - 4, . . . , 5K - 4, N} where 
5K — A < N < 5K + 1. 



Algorithm |3.1[ We estimate the three largest Lyapunov exponents Ai > A2 > A3 
and set A right = A 2 + 0.1(Ai - A 2 ) and A left = A 2 - 0.1 (A 2 - A 3 ). 



Algorithm |3.2 : As for Algorithm 3.1 



Algorithm 4.5: M = N, M' = 5, and d = (0, 1). 



Algorithm 4.6: Let M 1 = N, M[ = 5 and M 2 = N. 
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Figure 4: Comparing the approximations of the second Oseledets subspace W^fa) 
with the exact solution W 2 {x) which is known a priori. Each U N- approximation" is 
computed using cocycle data {A(T~ N x), A(T~ N+1 x), . . . ,A(T N x)}. The comparison is 



simply the Euclidean norm of the separation of the two unit vectors w 2 
and W2(x) G W 2 (x). 



(AO, 



x) e W. 



x 



Figure [4] compares the approximations yielded from the various approaches outlined 



in Sections^, ^ and [4] with the known solution of Equation (25). Each algorithm exhibits 



approximately exponential convergence with respect to the length of the sample cocycle 
up to (almost) machine accuracy of about 10~ 16 . Algorithm 4.3 is notably erratic whereas 
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the other algorithms converge smoothly, which suggests that in the limited data scenario 
(small N) it represents the less satisfactory choice. 



Algorithm 2.2 is slightly more accurate than the other algorithms for large N, while 



there is a limit to the accuracy of Algorithms 3.1 and 3.2 



It is worth noting that Algorithms 2^ and |4.6| do not perform as well when approxi- 
mating Oseledets subspaces corresponding to Lyapunov exponents A3, ... , A^. Whilst they 
reach machine accuracy with ease for W\ and W2, forming A(x,n) = A(T n ~ l x) ■ ■ ■ A(x) 
via numerical matrix multiplication produces greater inaccuracies for subspaces W3 
through Wi (see Remark [I]). On the other hand, Algorithms 3.1| 3.2 and 4.5 do not 
suffer from the same issue because they do not need to form A(x,n) but use only the 
generating matrices A(x). As such, they still reach machine accuracy, although a greater 
amount of data (larger N) is required. 
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Figure 5: Comparing the approximation of the second Oseledets vector with the exact 
solution for the two SVD based approaches, demonstrating the numerical instability which 
is overcome in the adapted SVD approach. 



Figure [5] is similar to Figure [4] except that it compares only Algorithms 2A_ and 2J2 



In doing so, it highlights the result of one of the numerical instabilities of Algorithm |2.1 
namely the pushing forward of [/■ (T~ N x) in Step 



Finally, Figure [6] shows the execution times of Algorithms |2.2[ |3.1[ |3.2[ |4.5| and |4.6 



which were timed using MATLAB's timing functionality. The most time-consuming step 



in Algorithm 4.5 is the SVD performed as part of the alterations from Section 4.1.1 



Algorithm 4J3 must perform two SVDs and Algorithm 2J2 must perform many more, 
which accounts for their longer execution times. 
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Figure 6: Comparing the execution time r of the various algorithms using MAT- 
LAB's timing functionality. Each algorithm is executed using the cocycle data 
[A(T~ N x), . . . ,A(T N x)}. 

5.2 Case Study 2: Particle dynamics - two disks in a quasi-one- 
dimensional box 

We consider the quasi-one-dimensional heat system studied extensively by Morriss et al. 
[2J [221 [271 Ell E2] which consists of two disks of diameter a = 1 in a rectangular box, 
[0, L x ] x [0, L y ], in which the shorter side, has length L y < 1o so that the disks may 
not change order. The two disks interact elastically with each other and the short walls, 
but periodic boundary conditions are enforced in the ^/-direction. The phase space of the 
system is then the set Icl 8 



X = (([0,Ljx[0,L,])/ 



~) 2 x 



X 



where ~ is the equivalence class associated with periodic boundary conditions, that is, 
(£2,2/2) e [0,^] x [0,L y ] have (xi,y{) ~ (^2,2/2) if Ui = IJ2 mod L y and x x = x 2 . 
The flow <f) T : X — > X consists of free-flight maps of time r, J- T : X — Y X, and collision 

maps C : X — > X so that <ft T (x) = C o J 77 "" o ■ ■ ■ o J 77 " 2 0C0 J 7 ^ (a;) where T\ H h r n = r. 

We consider a discrete-time version of the system by mapping from the instant after 
collision to the instant after the next collision, that is x i— > C o J^W(x) where t(x) is the 
free-flight time in the continuous system. The matrix cocycle is generated by the 8x8 
Jacobian matrices or the derivative of the flow evaluated instantly after each collision (i.e. 
A(x) = D (C o J 77 ^)) (x), see (5] for details). Due to a number of dynamic symmetries the 
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system has Lyapunov exponents Ai > A 2 > > — A 2 > — Ai with multiplicities 1, 1,4, 1 
and 1 respectively. This system has some symmetry, a high variation in expansion rates 
from iteration to iteration, a well separated spectrum, invertible Jacobian matrices, and 
relatively low dimension. Numerical integration of an orbit consisting of 4646 collisions 
yielded a sequence of Jacobian matrices {A(T~ 2323 x), . . . , A(T 2322 x)} which generate the 
cocycle A. 

For this model we use the same choice of parameters to execute the algorithms as 
with the previous model: 



Algorithm |2j2f M = N and {N k } = {1,6, 
5K — 4 < iV < 5K + 1. 



5k — 4, . . . , 5K — 4, N} where 



Algorithm 3.1 We estimate the three largest Lyapunov exponents Ai > A 2 > A3 
and set A right = A 2 + 0.1(Ai - A 2 ) and A left = A 2 - 0.1(A 2 - A3). 



Algorithm 3.2 : As for Algorithm 3.1 



Algorithm 4.5: M — N, M' = 5, and d = (0, 1). 



Algorithm 4.6: Let M 1 = N, M[ = 5 and M 2 = N. 



5.2.1 Criteria to assess the accuracy of estimated Oseledets spaces 

Since the Oseledets subspaces for this model are unknown, we test the approximations 
for two properties of Oseledets subspaces, namely their equivariance and the expansion 
rate, which defines the corresponding Lyapunov exponent. 



Equivariance: To test for equivariance, we approximate the second Oseledets 



vector, 



w 



(TV) 



(T n x) 



at each time 



n 



0,1, 



30. 



We then compute 



Af (a(x, n)w2 N \x) j — (T n x) and plot the result, where v&v/ \\v\\. If the ap- 



proximations are equivariant this value would be zero. 



Expansion Rate: To test the expansion rate, each approach is used to compute 

at time n = and we plot 



w. 



x 



W^ N \x), 



versus m. 



If W. 



(TV) 



X , 



is accurate, elements of W> 



should 



the second Oseledets vector, 
^log A(x } m)w { 2 N \x) 
grow at the correct rate: A 2 . 

Whilst the Oseledets vector u> 2 (x) must satisfy the above two properties, we must be 
careful when examining the results of these numerical experiments. For instance, (i) it is 
possible to choose vectors that are equivariant despite not being contained in any single 
Oseledets subspace, and (ii) any element of V 2 (x)\V 3 (x) 2 W 2 (x) (a much larger set than 
W 2 (x)) has Lyapunov exponent A 2 . 



5.2.2 Numerical Results 



Figure [7] shows the results of the equivariance test for the quasi-one-dimensional two disk 
model. At the lower end of cocycle data length (N = 75) all Algorithms except 3.2 



display reasonable equivariance, although Algorithm |3.1| remains equivariant for only a 
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handful of steps. For N = 150 and iV = 225 all approaches appear to produce close to 
equivariant results (note the changing scales in the vertical direction), with Algorithms 
3.1 and 3.2 lagging behind when N = 225. 
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Figure 7: The equivariance test for the various algorithms on the quasi-one-dimensional 
two disk model. Each approach is used to approximate the second Oseledets vector, 



w 



(T n x) G W 2 {N) {T n x) using cocycle data {A(T~ N x), ...,A 



[T x)}, at each time n 



0,1,..., 30. We then compute 



and plot the result. Note 



NA(x,n)w { 2 N \x) - w ( 2 N) (T n x) 

the different scales on each vertical axis. The plots shown are for N = 75 (top left), 
N = 150 (top right) and N = 225 (bottom). 



Figure [8] shows the results of the expansion rate test for the quasi-one-dimensional two 
disk model for various amounts of cocycle data {A(T~ N x), . . . , A{T N x)}. As expected, 
when there is a limited amount of data available (N small) the approximations either 
expand at the higher rate of Ai or only expand at the rate of A2 for a brief time before 
the error grows too large. As N is increased, the approximations expand at A2 for longer 
periods, suggesting that they more accurately represent w 2 (x). 

Most algorithms perform similarly regarding expansion rate. Note that the amount 
of cocycle data (size of N) needed to perform well in the Expansion Rate test is less than 
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Figure 8: The expansion rate test for the various approaches on the quasi-one- 
dimensional two disk model. The second Oseledets vector, ( x ) e 

(x), is approx- 
imated using cocycle data {A(T~ N x), . . . , A(T N x)} and we plot — log A(x, m)w < ^ r \x) 

versus m. If the approximation is accurate this quantity should tend to the value o 
A2 ~ 0.210, otherwise it would tend to the value of Ai ~ 0.325 both of which are shown 
in blue. The plots shown are for N = 25 (top left), A^ = 75 (top right) and A" = 150 
(bottom). 



that needed to perform well in the Equivariance test - this demonstrates the importance 
of good performance in both tests in order to assess whether or not the algorithms are 
performing well. 



5.3 Case Study 3: Time-dependent fluid flow in a cylinder; a 
transfer operator description 

An important emerging application for Oseledets subspaces is the detection of strange 
eigenmodes, persistent patterns, and coherent sets for aperiodic time-dependent fluid 
flows. In the periodic setting strange eigenmodes have been found as eigenfunctions of 
a Perron- Frobenius operator via classical Floquet theory; [25l |21| [26] . However, in the 
aperiodic time-dependent setting, Floquet theory cannot be applied. An extension to 
aperiodically driven flows was derived in [15], based on the new multiplicative ergodic 
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theory of [13]. Discrete approximations of a Perron- Frobenius cocycle representing the 
aperiodic flow are constructed and in this aperiodic setting the leading sub-dominant Os- 
eledets subspaces play the role of the leading sub-dominant eigenfunctions in the periodic 
forcing case. 

We review the four methods of approximating Oseledets subspaces with the aperiod- 
ically driven cylinder flow from [T3]. The flow domain is Y = [0, 2tt] x [0, tt], t G M + and 
the flow is defined by the following forced ODE: 

x = c — A(z(t)) sin(x — vz{t)) cos(y) + eG(g(x, y, z(t))) sin(z(t)/2) mod 2n 

y = A(z(t)) cos(x — vz{t)) sin(y). 

Here, z(t) = 6.6685zi(t), where zi(t) is generated by the standard Lorenz flow, 
A(z(t)) = 1 + 0.125 sm(y/5z(t)), G(ip) := l/(^ 2 + l) 2 and the parameter function 
if} = g(x, y, z(t)) := sxn(x—vz(t)) sin(y)+y/2— ir/4 vanishes at the level set of the stream- 
function of the unperturbed (e = 0) flow at instantaneous time t = 0, i.e., s(x, y, 0) = 7r/4, 
which divides the phase space in half. 




Figure 9: The second Oseledets subspace as determined by (a) Algorithm 2.2| (b) Algo- 
(c) Algorithm 3.2, (d) Algorithm 4.3 and (e) Algorithm 4.6 



rithm 3.1 



We set e = 1 as this value is sufficiently large to ensure no KAM tori remain in the 
jet regime, but sufficiently small to maintain islands originating from the nested periodic 
orbits around the elliptic points of the unperturbed system. 
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We construct the discretised Perron- Frobenius matrices P^\t) =: A(x) as described 



in Section 3 of [To], and briefly recapped in Example 1.2 , using a uniform grid of 120 x 60 
boxes, r = 8 and —32 < t < 32. In total, we generate 8 such matrices of dimension 
7200 x 7200. Thus, in this case study we have a limited amount of data, no symmetry, 
high dimension, and the matrices are non-invertible and sparse. 



In order to obtain reasonable results we executed Algorithms 2^ 3T 3.2 , 4.3 and 4J3 
with the following parameters: 



Algorithm 2.2: M = N = 4 and {N k } = {2,4}. 



Algorithm 3.1 We estimate the three largest Lyapunov exponents Ai > A2 > A3 
and set A right = A 2 + 0.1(Ai - A 2 ) and A left = A 2 - 0.1(A 2 - A 3 ). 



Algorithm |3.2 : As for Algorithm 3J. 



Algorithm 4.5: M' = M = 4 (so that only an SVD is used, and no push-forward 



step), N = 4 and c' = (0,1). 



Algorithm 4.6: M 1 = M[ = 4 and M 2 = 4. 



The results of these numerical experiments are shown in Figures [9] and 10 Recall that 
in this setting, the cocycle A(x,n) is a cocycle of discretised Perron-Frobenius operators 
acting on piecewise constant functions defined on Y; we identify these piecewise constant 
functions (with 7200 pieces) with vectors in IR 7200 . Figure [9] first shows the approximations 
of the second Oseledets vector w 2 (x) at time t = 0. In this setting the Oseledets vectors 



locate coherent structures: Figure 10 compares the push-forward of the approximations in 
Figure |9] with independently computed approximations of u; 2 (Tx) - the second Oseledets 
vector at time t = 8. 



In this study the data sample is insufficiently long for Algorithm 3.2 to work effectively. 



but the other algorithms produce similar results. A visual inspection of Figure 10 shows 
that the highlighted structures are approximately equivariant /coherent. 



6 Conclusion 



We introduced two new methods for computing Oseledets subspaces: one based on singu- 
lar value decompositions and the other based on dichotomy projectors. We also reviewed 
recent methods by Ginelli et al. [17] and Wolfe and Samelson [31], and presented an im- 
provement to both of these approaches that intelligently selected initial bases when only 
short time series were available to compute with. Finally, we carried out a comparative 
numerical investigation involving all four methods. 



Generally speaking, we found that Algorithms 2.2, 4.5, and 4.6 outperformed the 



dichotomy projector methods (Algorithms 3.1 and 3.2) when limited to moderate amounts 



of data were available, however, the dichotomy projector methods performed very well 
when long time series of matrices were available. The Ginelli approach (in particular the 



improved Algorithm 4.5) also worked very well with long time series. 



The improvements made to Algorithm 2.1 in Section 2.1 (namely the orthogonalisation 



step in Algorithm 2.2) produced an algorithm that could take advantage of longer matrix 
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Figure 10: Comparing the approximations of the second Oseledets vector w 2 (Tx) at time 
t — 8 with the push-forward of the approximations at time t = 0. Those labelled (a) are 
the push-forwards A(x, l)w^\x) whilst those labelled (b) are independently computed 
approximations w^\Tx) of W2{Tx). The algorithms used are as follows: (1) Algorithm 



2T2l (2) Algorithm 3Tl (3) Algorithm 3.2, (4) Algorithm 4.3 and (5) Algorithm KE 



sequences and return very accurate results. Of course, for each Algorithm one must 
choose the associated parameters sensibly to ensure good results. 

When only a short to moderate time series was available, we found mixed results in 
terms of the best algorithm. The improved SVD approach (Algorithm 2.2) was best for 
low to moderate length time series in the exact Toy model, while the improved Ginelli 
(Algorithm 4.5) and improved Wolfe (Algorithm 4.6) were marginally best in terms of 
equivariance and expansion rate, respectively for the 2-disk model. Each of these three 
algorithms produced similar results in the fluid-flow system. 

Choosing appropriate parameters for a particular application can be difficult. In the 
present review, good values were chosen by educated experimentation. On the other 
hand, the dichotomy projector methods, Algorithms 3.1 and 3.2, use parameters (A right 
and A left ) which can be chosen in a deterministic manner - by estimating Lyapunov 
exponents, which is a reasonably robust numerical procedure. Furthermore, a rigorous 
error approximation exists for Algorithm 3.2, a feature currently lacking for Algorithms 



2.2, 4.5 and 4.6 



The memory footprint of each approach scales quite differently with dimension. In 
Section 5.3, Algorithms 2.2 and 4.6 could take advantage of the sparseness of the d x 
d generating matrices of the cocycle. However, since A(x, n) is formed by matrix 
multiplication, for large n the matrix A(x, n) becomes dense and may require memory 
of the order of d 2 floating point numbers. The dichotomy projector Algorithms 3.1 and 



3.2, need to form an Nd x (N + l)d matrix, but with sparse generating matrices, this 
requires memory much less than of the order of d 2 floating point numbers. Algorithm 4.5 



has the most conservative memory footprint, but depends on its initialisation parameter 
M' and the Oselelets subspace number j. If M' is large, then A(T~ M x, M') in Step 1 
can become dense and require 0(d 2 ) floating point numbers. On the other hand, the 
stationary Lyapunov basis requires jd floating point numbers to be stored, so if j ~ d 
this can becomes comparable to d 2 . 

Section [5\3 involves non-invertible generating matrices and apart from Algorithm 3.2 



each approach succeeded in producing a reasonable solution, showing that the Algorithms 
can perform well in the non-invertible setting. Continuing with the non-invertible situa- 
tion, if one wishes to approximate Oseledets subspaces corresponding to negative numbers 
with very large magnitudes (Xj 



— oo 



then Algorithms 2.2 and 4.6 



may struggle as 

rapidly contracting directions (relative to the dominant direction corresponding to Ai) are 
quickly squashed during the matrix multiplication used to approximate A(x, n) leading 
to inaccurate numerical representation of A(x,n). 

The dichotomy projector approaches of Algorithms 3T and 3^ are able to compute Os- 
eledets subspaces corresponding to smaller, sub-dominant Lyapunov exponents A 3 , A 4 , . . . 
provided larger amounts of cocycle data is available. However, if Xj ~ — oo, we are forced 
to choose A nght or A left w — oo which means either problem (j8| or Q (in Algorithm 3.1 



which also feature in Algorithm 3.2) are ill-conditioned and fail. 

The same problem manifests itself in Algorithm 4J3 , even though it is able to compute 
Oseledets subspaces corresponding to smaller, sub-dominant Lyapunov exponents. The 
sum of the logarithm of the diagonal entries of the j x j generating matrices of the cocycle 
R(x, n) average to the logarithmic expansion rate of the j-parallelepiped formed at x by 



the stationary Lyapunov vectors s\ 



(oo) 



X) 



as it is pushed- forward. Thus, the 



logarithm of the zth diagonal entry of the generating matrices of R(x, n) has a time 
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average of Aj [I] and if \j ~ — oo, R(x,n) will feature diagonal entries close to, or equal 
to zero and R{x,n)^ 1 won't exist. 



In summary, Algorithms |2.2| and |4.6| are best suited to situations with limited cocycle 
data when one of the most dominant Oseledets subspaces is desired. Algorithm 4.5| 
can be applied to both limited and high data situations by choosing M' appropriately, 
and can compute most Oseledets subspaces provided their Lyapunov exponents are well- 
conditioned. If ample data is available and information regarding the system is lacking 
(making the choice of parameters for the other approaches difficult), the approaches of 
Algorithms 3A and 3J2 may be preferred for their relatively deterministic parameter 
selection. 
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